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Abstract. Numerical solutions of the incompressible magnetohydrodynamic (MHD) 
equations are reported for the interior of a rotating, perfectly-conducting, rigid 
spherical shell that is insulator-coated on the inside. A previously-reported spectral 
method is used which relies on a Galerkin expansion in Chandrasekhar-Kendall vector 
eigenfunctions of the curl. The new ingredient in this set of computations is the rigid 
rotation of the sphere. After a few purely hydrodynamic examples are sampled (spin 
down, Ekman pumping, inertial waves), attention is focused on selective decay and 
the MHD dynamo problem. In dynamo runs, prescribed mechanical forcing excites a 
persistent velocity field, usually turbulent at modest Reynolds numbers, which in turn 
amplifies a small seed magnetic field that is introduced. A wide variety of dynamo 
activity is observed, all at unit magnetic Prandtl number. The code lacks the resolution 
to probe high Reynolds numbers, but nevertheless interesting dynamo regimes turn out 
to be plentiful in those parts of parameter space in which the code is accurate. The 
key control parameters seem to be mechanical and magnetic Reynolds numbers, the 
Rossby and Ekman numbers (which in our computations are varied mostly by varying 
the rate of rotation of the sphere) and the amount of mechanical helicity injected. 
Magnetic energy levels and magnetic dipole behavior are exhibited which fluctuate 
strongly on a time scale of a few eddy turnover times. These seem to stabilize as the 
rotation rate is increased until the limit of the code resolution is reached. 
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1. Introduction 

In a previous paper, a spectral metliod for computing incompressible fluid and 
magnetohydrodynamic (MHD) behavior inside a sphere was introduced (Ref. [1], 
hereafter referred to as "MM"). The emphasis in MM was on accurate computation 
of global flow patterns throughout the full sphere (including the origin) with special 
attention to the computation of dynamo action, whereby the kinetic energy of 
a turbulent conducting fluid may give rise to macroscopic magnetic flelds. The 
computations were limited to moderate Reynolds numbers, with boundary conditions 
in which the normal components of the velocity fleld, magnetic fleld, vorticity, and 
electric current density were required to vanish at a rigid, spherical, insulator-lined, 
perfectly-conducting shell, and the three components of the velocity and magnetic fleld 
were regular at the origin. The feature previously lacking that we wish to explore in the 
present paper is that of uniform rotation of the sphere. We postpone to the future the 
investigation of insulating spherical shells which permit the magnetic fleld to penetrate 
into a non-conducting region outside [2] , though we remark later on some considerations 
relevant to this modiflcation. 

In Section [21 we formulate the equations to be solved for a uniform- density 
conducting fluid inside a sphere in a familiar set of dimensionless ("Alfvenic") units. 
We refer to MM for background and such of the details as remain unchanged. The 
main changes reported here are: (1) the introduction of a Coriolis term in the equation 
of motion (the centrifugal term may be absorbed in the pressure for incompressible 
flow); and (2) the velocity fleld v (instead of the vorticity lo) and magnetic fleld B 
are expanded in orthonormal Chandrasekhar-Kendall ("C-K") vector eigenfunctions of 
the curl [31 HJ [5], El [1]. The physical situation being simulated is again a perfectly 
conducting, mechanically impenetrable sphere coated on the inside with a thin layer of 
insulator, but now viewed from a coordinate system that is regarded as rigidly rotating 
with the bounding spherical shell. Unsurprisingly, the dynamical phenomena resulting 
are markedly different from, and richer than, they were in the case without rotation. 

Section [3] describes the results of some purely hydrodynamic runs (the code may 
be readily converted into a Navier-Stokes code by simply deleting the terms associated 
with the magnetic fleld). Included are examples [3, [H] of: (i) spin down, or decay of 
relatively rotating kinetic energy due to the action of viscosity; (ii) Ekman pumping 
with flow patterns that result from rotating boundaries; (iii) internal waves, three- 
dimensional relatives of meteorological Rossby waves, that depend on the stabilization 
introduced by rotation for their oscillatory features; and (iv) some mechanically-forced 
runs with a flnite angle between the symmetry axis of the forcing and the axis of rotation. 
This fourth case results in columnar vortices due to the effect of rotation [H [9] (note 
convection is absent in our present formulation). 

Section m proceeds to a consideration of the MHD case with an emphasis on selective 
decay and the kinds of dynamo behavior we have been able to resolve. The spectral 
method we use involves inherently less resolution than some other methods in use, and 
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we have been careful to study parameter regimes only where we can resolve the relevant 
length scales. Selective decay is observed to be somewhat arrested as the rotation rate 
is increased. A pleasant surprise has been the wide variety of dynamo behavior we have 
been able to resolve without the need to reach parameter regimes regarded as realistic for 
planetary dynamos [ini [HI El [121 [131 12] • In the summary, Section [3, we describe briefly 
some plans we have for improving the resolution of the code by some pseudo-spectral 
modifications and some intended future diversification of the boundary conditions. 

2. Computational method 

We begin from the magnetohydrodynamic (MHD) equation of motion in a rotating 
coordinate frame p3l [2] , 



In the dimensionless Alfvenic units [T], v is the vector velocity field, B is the magnetic 
induction, a; = V x v is the vorticity, and j = V x B is the electric current density. 
The generalized pressure is V. The dimensionless viscosity, which, in the dimensionless 
variables, can be interpreted as the reciprocal of a mechanical Reynolds number, is u, 
and the magnetic diffusivity, which can be interpreted as the reciprocal of a magnetic 
Reynolds number, is rj. The vector field f is a solenoidal, externally-applied forcing 
field which is intended to mimic the presence of mechanical sources of excitation 
of V. Equations ([T|) and ([2|) are to be supplemented by the requirements that the 
divergences of both v and B must vanish everywhere. Dropping equation ([2]) and the 
terms in equation ([T|) containing j and B leaves the forced Navier-Stokes equation, and 
dropping f leaves the unforced version of Navier-Stokes. Q is the (constant) rotation 
speed of the coordinate system, understood to be attached to a rotating spherical shell 
that constitutes the boundary and is both mechanically impenetrable and perfectly 
conducting, with a thin layer of insulator on the inside surface. 

The non-trivial boundary conditions imposed are that the normal components of 
V, B, j, and a; shall all vanish at the radius of a unit sphere centered at the origin. The 
three components of the fields v and B are also required to be regular at the origin. 
The vanishing of the normal components of v and uj at the surface of the unit sphere 
are implied by, but do not imply, no-slip boundary conditions at that radius. Going 
further with an attempt to implement fully a set of no-slip boundary conditions raises 
unresolved paradoxes with respect to the pressure determination which we prefer not 
to confront here (see Refs. [15l [161 [I] for a discussion of these), believing that their 
seriousness and intractability require consideration in the context of simpler situations 
than the present one. 




(1) 



and the MHD induction equation. 




(2) 
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The spectral technique implemented involves expanding v and B in terms of C-K 
functions (defined below): 



v(r,t) = EQ™WJ./-W, (3) 

and 



■aqlri 

qlm 



Bir, t)=J2Qm{t) J qimir). (4) 

qlm 

The C-K functions P H El El [1] J i are defined by 

J. = AV X r^i + V X (V X ripi) , (5) 

where we work with a set of spherical orthonormal unit vectors (f, 6, (p) and the scalar 
function ipi is a solution of the Helmholtz equation, (V^ + X'^)ipi = 0. The explicit form 
of ipi is 

i:,{r,9,^) = C,iji{\Xgi\r)YiU0A), (6) 

where ji{\Xqi\r) is a spherical Bessel function of the first kind which vanishes at r = 1 
and Yim{0,(f)) is a spherical harmonic in the polar angle 6 and the azimuthal angle (p. 
The subindex i is a shorthand notation for the three indices {q,l,m), where q indexes 
the successive values of A that make ji vanish at r = 1 for each value of /; g = 1, 2, 3, . . . 
corresponds to the positive values of A, and q = —1, —2, —3, . . . indexes the negative 
values; finally Z = 1, 2, 3, . . . and m runs in integer steps from — / to +/. The vectors Jj 
satisfy 

V X Ji = AiJi, (7) 

and with the proper normalization constants are an orthonormal set that has been shown 
to be complete [6]. The integral relation expressing the orthogonality of the Jj is: 

Jg/m. ■ Jg',/',™'*^^ = Sq^qi6l^li6m,m', (8) 

where the asterisk denotes complex conjugate, and with the normalization constants 
given by: 

Cqi = \Xqiji^^{\Xqi\)\-' [l{l + (9) 

The scheme for solving equations ([I]) and ([2]) is conceptually simple. We substitute 
the expansions ([3]) and (jl]) into equations ([1]) and ([2]), utilize the fact that the Jj are 
eigenf unctions of the curl, and then take inner products one at a time with the individual 
Jj. Their orthogonality enables us to pick off expressions for the time derivatives of the 
time dependent expansion coefficients ^'^ and ^f, and equations ([T]) and ([2]) are thereby 
converted into a set of ordinary differential equations for the expansion coefficients. 
These appear as 

^ = E 4fc - ) +2j2n.o]- u)^a + e/, (10) 
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Figure 1. (a) Time liistorics of the decay of meclianical energy for four liydrodynamic 
runs (HI to H4) witli identical initial velocity fields but different rotation rates fi. (b) 
Decay rates vs. the square root of for the same four runs. 



and 

^ = E55.eKf-^A,^ef, (11) 

with the couphng coefficients defined as 

= ^kl]k^ Bj,^ = (12) 

p.^ = J J*. 3 J X 3kdV, O) = J3*x 3,dV. (13) 

The infinite set of ordinary differential equations is truncated at some level above 
maximum values of |g| and /, in the usual manner of a Galerkin approximation |17j . 
The evaluation of equations (llUp and (II ip and the storage of the resulting arrays of 
coupling coefficients in tables, are the most demanding numerical tasks of the problem. 
Once available, they do not have to be recomputed, and provide a method for verifying 
the ideal quadratic conservation laws with high accuracy [Tj. Also, since the method 
is purely spectral and fields are only computed in real space for visualization purposes, 
there is no numerical singularity in the center of the sphere. 

The main drawback of the scheme, as with any wholly spectral one, is that the 
convolution sums in equations (fTOl) and ( iTTl) grow rapidly with increasing maximum 
values of |g| and /, and limit the resolution when compared to pseudospectral 
computations utilizing fast transforms (in practice, a resolution of max{|g|} = max{/} = 
9 was used in all the runs). This limits us to modest Reynolds numbers (all our 
computations reported here have limited themselves to resolvable Reynolds numbers). 
Future plans include pseudospectral modifications to the evaluation of at least the 
angular parts of the nonlinear terms in equations ([1]) and ([2]), as will be mentioned 
again in the final summary (Section [5]). 
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Figure 2. (a) Enstrophy for runs HI to H4 as a function of time (same labels as in 
Figure [T|). (b) Energy spectra at t = 1 for the same runs, and (c) energy spectra at 
t = 6. In all cases, solid lines are for il = 0, dotted lines for fl = 1, dashed lines for 
= 4, and dash-dotted lines for ft = 10. 

3. Hydrodynamic examples 

Some neutral-fluid effects (good introductions to all of which may be found in Refs. 
[TllH]) are treated first before proceeding to MHD. It is worth noticing here that although 
our boundary conditions are implied by, but do not imply no-slip velocities, several 
qualitative and some quantitative agreements are observed with previous experiments 
and theory. First we study simple problems in which initially relatively rotating fluids 
adjust themselves to rigid rotation with the spherical shell. We study these decays as 
functions of the rotation rate Q. The prediction [7] is that the decay of the non-rigid body 
components should be exponential, with a decay rate that varies as gure [T]^a) 

shows the time histories of the decay of mechanical energy for four runs with identical 
initial velocity fields limited to a few random low mode numbers (large spatial scales). 
The runs are delineated as runs HI, H2, H3, and H4. Specifically, we have as initial 
conditions a random superposition of modes with |g| = 1,2, / = 1,2, and all allowed 
values of m, with viscosity u = 0.01 and values of Q of 0, 1, 4, and 10 respectively. Each 
curve is approximately exponential, and when they are fitted with exponentials, the 
decay rates plotted vs. the square root of fl appear as in Figure [T](b), and are adjudged 
to be in satisfactory agreement with theory [7j. The nonlinearities excite smaller spatial 
scales, and the decay process is progressively enhanced by increasing the rotation rate. 

A second feature of rotating spherical fiows observed in our code is the development 
of Ekman-like layers and the action of Ekman pumping [71 [8] (see also [18] for a 
detailed study in rotating spherical shells). The fiow patterns are characterized by 
the development of interior vortical fiows with some symmetries, and thin layers that 
separate the large vortices and also lie along the wall boundary layers. These have a 

1/2 

characteristic thickness of the order of 5 ~ Ej^ R, where i? = 1 is the radius of the 
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Figure 3. Above: Mechanical energy density and velocity field lines in run El, at 
< = (left) and at t = 6 (right). Below: idem for run E2. For convenience, energies 
and field lines are always shown in pairs, with energy densities on the left and field 
lines on the right. The field lines change color according to the distance integrated 
from the initial point, from red to yellow, blue, and magenta. The red, green, and 
blue arrows indicate respectively the x, y, and z axis, is in the z direction. In both 
cases, the energy density is symmetric with respect to the equator, while the flow itself 
is antisymmetric in run El (above). 

sphere, and the Ekman number is Ex = vVL~^L~^, with L a characteristic scale of 
the flow. The ability of the code to compute these layers is limited by its resolution. 
Realistic values of the Ekman number are, for planetary core regimes [TH |9l [2], beyond 
our range. In all the runs presented here we will limit ourselves to cases where 5 can be 
properly resolved with the number of modes used in the simulations. The presence of the 
Ekman layers is concomitant with the development of smaller spatial scales and hence 
of more rapid dissipation. Figure [2] illustrates this fact. Figure EJj^a) is the integrated 
squared vorticity, or enstrophy, for the four runs whose decay has just been seen to be 
exponential. The largest rotation rate corresponds to the curve with the highest early 
peak in enstrophy spectrum, the second highest with the second largest, and so on. 
Figure [21(b) shows the energy spectra at t = 1 for the four runs and Figure Wis) shows 
the same energy spectra at t = 6. In every case, the flatter spectra, and hence the 
shorter wavelength dominances, correspond to the higher values of f2. This is the result 
of the formation of a thinner boundary layer as VL is increased. 

The range of Ekman layer behavior we have been able to observe is very wide. We 
show in FigureOthe results of two simulations (labeled El and E2) with initially random 
axisymmetric (m = 0) velocity fields which are purely azimuthal. For both runs, the 
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Figure 4. Equatorial cross section of the velocity field for an inertial wave in the 
rotating sphere. The reference frame is fixed to the sphere. H is in the z direction. 
The axial velocity is indicated by the colors, while the radial and azimuthal velocities 
are indicated by the arrows. For this mode, the frequency is w f« 49. 



initial velocity is proportional to the difference between Jg,i,o and J-q^i^, which has only 
an azimuthal component. For run El, q = 1 and / = 1, and for run E2, q = 1 and / = 2. 
Both runs have u = 0.01 and Q = 10. The time evolution of the runs is similar to the 
evolution displayed in Figures [1] and [21 However, the axisymmetric initial conditions in 
runs El and E2 make visualization of flow patters easier. Figure [3] shows the initial and 
late-time flow patterns for these runs, using the VAPOR graphics package [19] that will 
be repeatedly used throughout this paper for graphical demonstrations. The rotation 
generates poloidal components of the velocity field fast, and at late times different 
patterns are observed depending on the initial value of I. In run El, at late times the 
flow displays a poloidal circulation on top of the initial toroidal fleld: the flow is directed 
towards the center of the sphere along the axis of rotation, and a return flow is observed 
in both hemispheres close to the wall. In other words, the flow can be described as the 
superposition of a toroidal differential rotation and a poloidal meridional circulation. 
This circulation is radially outward in the meridional plane, directed towards the poles 
close to the wall, and redirected toward the equatorial plane again as the flow gets close 
to the poles. Both hemispheres show the same pattern. In run E2 the pattern is more 
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Figure 5. (a) Magnetic energy Em (solid line), kinetic energy Ey (dotted line), and 
magnetic lielicity Hm (dashed line) in the selective decay run SI. (b) Same quantities 
for run S2. The two power laws are indicated in the figures only as a reference. 

complex, and vertical velocities are observed in the vicinity of the axis of rotation, while 
high and at intermediate latitudes a poloidal circulation is generated. 

As a third hydrodynamic test of the code we demonstrate inertial wave motion 
in the small- amplitude limit. Equations ( ITOl) and ( fTTl) can be linearized in powers 
of a small departure from a uniform rotation velocity; then solutions can be sought 
which vary with time as e*™*. The resulting linear homogeneous algebraic system can 
be solved in a Galerkin approximation by expanding the velocity and vorticity in the 
C-K functions. An anti-Hermitian matrix results whose eigenfunctions can be found 
numerically and whose corresponding eigenvalues iw may be computed numerically in 
the process. Then any one of the oscillatory modes can be loaded numerically into the 
ideal version of the code and run with the overall amplitudes chosen to be very small. 
The time evolution is accurately predicted by the computation of the single modes, 
which are standing waves. Figure H] shows four equatorial cross-sections of the sphere at 
different times with the axial velocity indicated by the color codes, and the radial and 
azimuthal velocities indicated by arrows {Q = 10 in this run). The times are chosen to 
be one quarter-period apart. The oscillation frequency w for the particular mode shown 
as obtained from the eigenvalue problem is in good agreement with the results obtained 
from the fully non-linear code [w ~ 49). 

In forced hydrodynamic simulations in the presence of strong rotation, we also 
observe the development of columnar structures in the flow, aligned with the axis of 
rotation. The discussion of these simulations will be left for the next section, where the 
connection between these columns and dynamo action will be considered. 
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Figure 6. Relative helicity Hm /Em in selective decay runs SI (solid) and S2 (dashed). 



4. MHD and the dynamo 

Selective decay in the sphere 

Before passing to a discussion of the mechanicaUy driven spherical dynamo, we present 
first the resuhs of two tests of three-dimensional MHD "selective decay" as affected by 
the presence of rotation. Selective decay is a familiar turbulent decay process, usually 
incompressible, long studied in periodic geometry [2ni[211[22], wherein one ideal invariant 
is cascaded to short wavelengths and dissipated while another remains locked into long 
wavelengths and is approximately conserved. The phenomenon, closely connected with 
inverse cascade processes for driven systems, leads toward a state in which the ratio 
of the two ideal invariants involved is minimized and which therefore is accessible to 
variational methods. In 3D MHD, a quantity that may be preferentially dissipated 
is the total energy E = Ey + Em (kinetic plus magnetic) while magnetic helicity 
Hm may be approximately conserved. Under other circumstances, energy may be 
dissipated while cross helicity K is approximately conserved, leading to the phenomenon 
of "dynamic alignment," [231 l2ll [25] in which the velocity field and magnetic fields are 
highly correlated. The definitions of Hm and K are 



where A is the vector potential whose curl is B, and the integrals run over the entire 
volume of the fluid. What we are interested in demonstrating here is the effect that 
rotation has on the development of the selective decay of total energy relative to 
magnetic helicity inside a sphere. It will be useful for these and other purposes, to 
have a definition of the magnetic dipole moment: 




(14) 




(15) 




(16) 
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Figure 7. (a) Trace of the dipolc moment on the sm-facc of the unit sphere (above) 
and amplitude of the dipole moment as a function of time (below) for the selective 
decay run SI. (b) Same quantities for run S2. 

which is seen to be readily expressible in terms of the expansion coefficients for B [1] . 

The initially excited modes for the two runs we will present ("SI" and "S2") are 
those for q = ±3, / = 3, and all possible values of m. The initial values chosen for the 
expansion coefficients are: 

^±3,3,0 = -^0, ^±3,3,0<m<3 = ^o(l + i) , (17) 
^3,3,0 ~ "^^-3,3,0 ~ ^05 ^3,3,0<m<3 ~ "^^-3,3,0<m<3 ~ ^o(l ^ (l^) 

with uq and Bq chosen so that at t = 0, the magnetic and kinetic energies are 
Em = Ey ^ 0.5, K = 0, and Hm ~ 0.034. Some helicity cancellation occurs because 
of the two signs of A (or q). As a comparison, note that for the g = 3, / = 3 mode 
alone, Hm/Em is no more than about 0.072 (this is the maximum value of \Hhj/Em\ 
if only modes with |g| = 3, / = 3, and one sign of A are excited). In both runs, the 
magnetic diffusivity and kinematic viscosity are v = rj = 0.006; the Reynolds numbers 
are Re ~ Rm ~ 170, based on the radius of the sphere. The two runs differ by the values 
of Vt chosen, which are 2 and 12, respectively. These mean that the Rossby and Ekman 
numbers of the two runs are, respectively, Rq = U{QR)~^ = 0.5, E^ = 0.003 for SI and 
Ro = 0.083, Ek = 0.0005 for S2. The decay of magnetic energy, kinetic energy, and 
magnetic helicity for Runs SI and S2 are shown in Figure [51 The behavior in Run SI 
is not significantly different from the non-rotating case [1]. Note that in both runs, the 
kinetic energy at late times is negligible, and that magnetic and kinetic energies decay 
faster than the magnetic helicity. However, the decay of all these quantities in S2 seems 
to be faster than in SI. 

The relative helicity, Hm/Em, is shown for Runs SI and S2 in Figure [61 It will be 
seen that the increased rate of rotation in S2 has somewhat arrested the selective decay. 
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Figure 8. Above: kinetic energy density and velocity field lines (left) and magnetic 
energy density and magnetic field lines (right) in the initial conditions of selective 
decay runs SI and S2. Middle: same fields at t = 14 in run SI. Below: same fields at 
t = 14 in run S2. Colors and labels are as in Figure O 

for reasons not totally understood. It may be that the rotation has resulted in sufficient 
two-dimensionalization of the flow [TJ El E] that the inherently three dimensional nature 
of the selective decay has been compromised. But from Figure [6] it can also be seen 
that the ratio of Hm to Em has approached reasonably closely to its maximal value of 
min~^{|Agi|} ~ 0.22 (the maximum value of \Hm/E]\i\ when only modes with |g| = 1, 
/ = 1, and one sign of A are excited). The maximal value would indicate a total 
disappearance of the non-rotational kinetic energy (i.e., a rigid rotation), and all the 
magnetic energy in the largest-scale modes (smallest |A|) allowed by the boundary 
conditions: a magnetized, "frozen" condition. 

The behavior of the dipole moment for the two runs is shown in Figure [3 The solid 
lines are the traces of the projected direction of the dipole moment on the surface of 
the sphere as functions of the time, and the orientation is such that the axis of rotation 
points upward. In the lower parts of these two figures, the magnitude of the dipole 
moment is plotted vs. time. In both cases, it will be seen that the dipole's orientation 
initially wanders erratically near its initial position, and finally ends at a "mid-latitude" 
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Table 1. Dynamo runs: q and I give the scales where mechanical energy is injected by 
the forcing, is the rotation rate, and v and rj are respectively the kinematic viscosity 
and magnetic difTusivity. "Helical" indicates whether the forcing injects mechanical 
helicity, and "Axisym." indicates whether the forcing is axisymmctric. a is the angle 
between Jl and the z-axis (the axis of symmetry in axisymmctric forcings). Finally, 
Ro, Ek, and Re are respectively the Rossby, Ekman, and Reynolds numbers, all based 
on the radius of the sphere. A resolution of max{|(7|} = max{Z} = 9 was used in all 
the runs. 
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direction not far from where it began. This was something of a surprise to us, since 
we had expected it to hne up with the axis of rotation or at least close to it. Figure [8] 
shows VAPOR plots with the energy densities, velocity and magnetic field line structure 
for the initial conditions in Runs SI and S2, as well as the late stages of both runs 
[t = 14) when the selective decay process has saturated and all the nonlinear terms are 
small, preventing any further evolution of the system except for dissipation. For strong 
rotation (run S2), the velocity field is quasi- two dimensional and develops column-like 
structures at late times, while the magnetic field is highly anisotropic (although the 
dipole moment is not aligned with the axis of rotation). Note magnetic field lines in 
this case are aligned with the z axis (the axis of rotation), and velocity field lines are 
mostly toroidal. Actually, the ratio \v^/vz\ at t = 14 averaged over the whole volume 
for this run at t = 14 is ~ 13. 

4-2. Dynamos 

We turn now to the case of the forced spherical dynamo computations, in which specified 
non-zero forcing functions f are added to the right hand side of equation ([T]) or (fTOll to 
provide a persistently-active, non-decaying velocity field. After a purely hydrodynamic 
run to reach a statistically steady state, very small magnetic fields are introduced to 
see if the velocity fields will cause them to amplify, and attention focuses on questions 
like the orientation of the resulting magnetic dipole moment relative to the axis of 
rotation and the dipole moment's magnitude. We are also interested in the kinetic and 
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Figure 9. Left: kinetic energy density and velocity field lines at late times in run Df , 
when the dynamo has saturated. Right: magnetic energy density and magnetic field 
lines at the same time. Colors and labels are as in Figure [3l 



magnetic energy spectra that result, and how the eponymous dimensionless numbers of 
the rotating fluid (Reynolds, magnetic Reynolds, Rossby, Ekman) influence the magnetic 
quantities. These appear to be essential control parameters of the problem. The range 
of possibilities is clearly very wide, and we have not begun to explore the entire space 
of possible parameters. Rather, we content ourselves with showing samples of different 
behavior that have emerged for different combinations that lead to regimes which the 
code will resolve satisfactorily. Our efforts should not be compared with explorations 
of the space of parameters in realistic geodynamo simulations (see e.g. [26]) but rather 
as an extension of dynamo simulations of incompressible MHD flows (often done using 
periodic boundary conditions [27l [28l [29l [30]) to include the effect of boundaries and 
rotation. 

Several of the forcing functions used in the runs we will display are axisymmetric, 
but their axes of symmetry are not aligned with the axis of rotation. The resulting 
overall asymmetry quickly excites all the available retained modes, to some degree. The 
general form of the axisymmetric forcing function used is 



For any value of A and B, this superposition of C-K functions gives an axisymmetric 
forcing f. For A = B = 1, there will be no net helicity involved in the forcing (the 
curl of f is perpendicular to f), and the only non- vanishing component of velocity 
that is forced is the (j) component (the azimuthal component with respect to the axis 
of symmetry of the forcing, the z axis). This forcing can be considered as a simple 
differential rotation, where the number of nodes in v^{r, 9) is controlled by the values 
of q and /. The axis of rotation is typically oriented at some specified angle a (often 
30°) to the forcing function's axis of symmetry (the polar axis, in spherical coordinates). 
Thus the rotational motion and the forcing have no shared symmetry, and the resulting 
mechanical motion is totally asymmetrical. 

Non-axisymmetric forcing functions are obtained by superposing C-K modes with 




(19) 



ql 
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Figure 10. Left: kinetic (dashed black line) and magnetic (solid blue line) energy as 
a function of time in run Dl. Right: cosine of the angle between H and /x in the same 
run. 



m ^ 0, i.e. 

f = Xl^gimJg^m, (20) 

qlm 

and when ^^^^^ = C,^^ i ^ the forcing is non-hehcaL In any other case the curl of f has a 
projection into f , and the forcing injects mechanical helicity into the flow. 

In an effort to systematize the runs we have done and the reasons we have done 
them, we have assembled the important parameters for the dynamo runs (labeled Dl 
through D12) in Table [TJ Listed in Table [1] are: the q and / values where the forcing 
was concentrated (determining its characteristic length scale); the rigid rotation rate; 
the kinematic viscosity (reciprocal Reynolds number, if the kinetic energy is close to 
unity) and magnetic diffusivity (this study is restricted to the magnetic Prandtl number 
Pm = ly/rj = 1 case); an indication of whether the forcing was axisymmetric or not; the 
angle between the axis of rotation and the axis of symmetry of the forcing, when the 
forcing function has an internal symmetry; the Rossby and Ekman numbers of the flow 
into which the seed magnetic field is introduced; and whether or not the forcing injected 
net mechanical helicity. 

It is perhaps worthwhile to say a word about the motivation for the progression 
of runs shown in Table [H The first remark is that it seems to be relatively easy to 
excite a dynamo and generate a dipole moment, but relatively difficult to generate one 
that behaves according to our predispositions and hopes: a dipole moment with some 
alignment with the axis of rotation, and that reverses periodically or randomly with 
long times between reversals (compared with the turbulent turnover time). We have 
found wild oscillations in both magnitude and direction that seem to decrease with 
decreasing Rossby and Ekman numbers. Since the Reynolds numbers are limited by the 
resolution, the principal means of decreasing both the Rossby and Ekman numbers is by 
increasing the rotation rate Q. This, however, eventually decreases the thickness of the 
boundary layer below the resolution of the code, and beyond that point, the accuracy 
of the computations becomes suspect. The lowest Rossby and Ekman numbers that 
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Figure 11. (a) Magnetic energy in dynamo runs Dl (solid), D2 (dotted), D3 (dashed), 
D4 (dash-dotted) and D5 (dash-triple dotted), (b) Magnetic energy spectrum for the 
same runs at late times, after non-linear saturation of the dynamo takes place. The 
arrows on top indicate the scale where mechanical energy is injected in each run; from 
left to right: Dl, D2, and D3-D5. The magnetic energy spectrum corresponding to 
Run Dl has been multiplied by a factor of 100. 



appear in the entries of Table [T] represent those below which the resolution limitations 
are encountered. It will be seen below that as we progress toward them, the dipolar 
behavior looks gradually more like what we might expect, the dipole moment gets 
stronger and more aligned with the axis of rotation, and the time between reversals gets 
larger. 

We begin by showing some results for a weak dynamo situation, Run Dl, where 
the magnetic energy Em remains always much smaller than the kinetic energy Ey- In 
the forcing function (flQll . the driven modes have q = l = l, A = B = l,i' = ri = 0.002, 
a = 30°, and fl = 2. The amplitude of the forced modes is 1^^ / qI — O-^- The two 
Reynolds numbers, and Rm, based on the measured r.m.s. velocity before the 
magnetic seed is introduced and the radius of the sphere, are both about 500. The 
Rossby number is Rq = 0.5, and the Ekman number also based on the radius of the 
sphere is E^ = 1 x 10"^. Figure [9] shows the streamlines of the flow, magnetic field 
lines, and the energy densities at late times, once the dipole is established {t ~ 80). In 
the steady state, the magnetic dipole moment |^| is of the order of 0.001 and the ratio 
of magnetic to kinetic energies is Em/Ev ~ 0.0005. The magnetic energy rises to a 
characteristic value and oscillates somewhat irregularly as shown in Figure [TOl while the 
cosine of 7 (the angle between the axis of rotation and the orientation of ^) oscillates 
with roughly the same periodicity (see also Figure [TOil with almost a 180° variation in 
some cases, and with an average orientation almost perpendicular to Q. In this weak 
case Dl, the hydrodynamic flow remains laminar, stable, and almost time independent. 

The global evolution of the system is similar to what we will show in the remaining 
runs. Once the magnetic field is introduced at t = 0, and if Rm is large enough, 
the magnetic field is amplified exponentially (this stage is often called the "kinematic 
dynamo" regime) until the Lorentz force modifies the flow and non-linear saturation is 
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Figure 12. Left: (a) Trace of the dipole moment on the surface of the umt sphere 
(above), ampUtude of the dipole moment (middle), and cosine of the angle between 
the dipole moment and the axis of rotation as a function of time, for run D3. (b) Same 
quantities for run D5. (c) Same quantities for run D7. 



reached. At late times, an MHD state is reached in which magnetic energy is sustained 
against Ohmic dissipation by dynamo action. In run Dl, the flow is reminiscent of 
the hydrodynamic flow previously described in Figure [31 Although the forcing is 
axisymmetric and purely toroidal, rotation generates a poloidal circulation and as a 
result the flow points outwards in the equatorial plane, and inwards along the axis 
of rotation. Each hemisphere has mechanical helicity of opposite signs, while the net 
mechanical helicity of the system fluctuates around zero. The magnetic field seems to 
be sustained by an a-Q mechanism, where the differential rotation is sustained by the 
mechanical forcing and the a-effect is given by the Ekman-like circulation. Magnetic 
energy is concentrated in the center of the sphere, where the flow has a stagnation point. 

A considerably stronger dynamo than Dl is represented in Run D2, where in the 
excitation function f|T9|) we choose q = I = 2 and l^g^ol ~ Again, A = B = 1, 
v = 7] = 0.002, and f2 = 2. The magnetic moment rises from zero, and attains a typical 
magnitude of ^ 0.5, about 100 times larger than in Dl. The Reynolds numbers are 
Re = Rm ~ 500, and the ratio of magnetic to kinetic energy oscillates around 0.15. The 
cosine of 7 (the angle between fi and ft) oscillates wildly in time, and the orientation 
of the dipole shows no preferred direction. The only difference between run Dl and D2 
is the change in the forcing scale, and the result seems to indicate a separation of scales 
between the forcing and the largest scale in the system helps the dynamo, as indicated 
by the larger ratio Em/Ev in run D2, and as also reported before in simulations with 
periodic boundary conditions [28], [31], [32] . In all these runs, the largest available scale is 
fixed and given by the inverse of the smallest |A| (corresponding to A±i^i) and determined 
by the radius of the sphere {R = 1), while the separation between this scale and the 
forcing scale is controlled by the values of q and / in the forcing function (see Table [T|). 
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Figure 13. Magnetic energy in dynamo runs Dl (solid), D2 (dotted), D3 (dashed). D5 
(dash-dotted) and D7 (dash-triple dotted). Note the intermittent growth of magnetic 
energy at early times in Run D7. 

The largest the values of |g| and /, the smallest the scale where mechanical energy is 
injected. 

Several runs were done (see e.g. the runs Dl to D5 in Table [1]) in which the forcing 
was gradually moved to smaller scales (from g = / = linDltog = / = 3in D3), and 
in which the rotation rate was progressively increased (from Q = 2 in runs D1-D3 to 
r2 = 8 in D5). As these changes were made, the amplitude of the forced modes |^^; ^| in 
equation ( fT9l) had to be increased in order to reach a statistically steady state with r.m.s. 
velocities of order one before the magnetic field was introduced (their amplitudes were 
0.4, 1.1, 1.6, 2.2, and 3.6, from run Dl to D5). The reason for this can be understood 
as follows: the Coriolis force in equation ([T]) acts as a restoring force that opposes the 
growth of perturbations. This is also the reason why this system can sustain waves, as 
was shown in Section [3l In all these runs. A, B, u, rj, and the angle of inclination a 
were kept the same. As will be seen from Table 1, all the forcing was non- helical. 

All five runs are considered to have been able to resolve the Ekman layers that 
developed, but they would likely not have been resolved at higher values of Q. Figure [TT] 
shows the general trend resulting from the smaller scale forcing and increased rotation. 
Figure fTlT a) shows a logarithmic-linear plot of the total magnetic energy versus time for 
runs Dl to D5. The five runs showed increasingly large growth rate, a higher saturation 
level of Em/Ev, and increasing In Run D5, the ultimate ratio of Em to Ey was 
about 0.4 and \iJ,\ was close to unity. Figure fTTT b) shows magnetic energy spectra for 
runs Dl to D5, with decreasing Rossby number Rq (based on the radius of the sphere) 
of 0.5 (runs Dl to D3), 0.25 (D4), and 0.125 (D5). It will be seen that there develops a 
small excess of magnetic energy in scales larger than the forcing scale with decreasing 
Rossby numbers. Note the development of a "bump" in the magnetic energy spectrum 
at A ~ 9 in runs D4 and D5. 
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Figure 14. Above: kinetic energy density and velocity field lines at late times in Run 
D5 (left), and magnetic energy density and magnetic field files at the same time for the 
same run. Below: same quantities for Run D7. Note the development of anisotropics 
in the presence of large fl in this run. 

A second trend is indicated in Figure [T2j namely the orientation of the dipole 
moment onto the axis of rotation seems less erratic with decreasing Rossby and Ekman 
numbers. That is, there are more eddy turnover times (in units oi R/U) between the 
near reversals as Rq and Ek are decreased, and the projection of the dipole moment 
on the unit sphere gets more localized around the two "poles" defined by the axis of 
rotation. We borrow here the term "reversal" from the palaeomagnetic record, where 
during a reversal the orientation of the dipole moment changes about 180° and its 
amplitude decreases, in opposition to an "excursion" in which the direction and the 
orientation of the dipole moment changes in a short period of time without resulting in 
a full reversal [33] . 

To verify this behavior, in Runs D6, D7, and D8 we successively decreased Rq and 
Ek while keeping the other parameters constant. At the present resolution, we could 
not decrease the values of Rq and E^ below the values for run D8 while keeping the 
boundary layer well resolved. In geodynamo simulations, a similar effect was reported, 
and it was noted that the behavior of the dipole moment was controlled by the amplitude 
of the Rossby number, independently of the values of the Ekman and Rayleigh numbers 
[26]. 

Figure [12] shows the trace of the dipole moment fi on the surface of the unit sphere, 
its amplitude, and the cosine of the angle between fj, and ft for runs D3, D5, and D7. 
While in Run D3 the trace of fx fills the entire surface of the sphere, as Rq and E^ 
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Figure 15. Kinetic energy density (left), velocity field lines (middle), and view from 
top of the kinetic energy density superposed with velocity field lines (right) in nm DIO, 
before magnetic energy is introduced. Note the columnar structures in the velocity field 
aligned along n (in the z direction). 

are decreased seems to fluctuate around two regions in opposite sides of tlie spliere. 
Tliese regions get more localized with decreasing Rq and Ek. Also, the time between 
excursions of /j, outside these regions gets larger, as shown by cos(7). In Run D7, after a 
transient that finishes at t ^ 80, cos(7) stays at 1 or —1 for ^ 20 turnover times before 
changing sign rapidly. It is also worth noticing that as Rq and decrease, the time it 
takes the system to develop a dipole moment of order one gets larger (see the evolution 
of at early times in Figure [T2l) . This is also observed in the time evolution of the 
magnetic energy (see Figure [T3l) . Instead of having an exponential growth of Em at early 
times as in Runs Dl to D5, Runs D6 to D8 show a more intermittent behavior: magnetic 
energy grows in rapid "bursts" and stays around that value until a new burst increases 
the magnetic energy again. In Run D7, this process saturates around t ~ 80 and no 
further change in the average value of Em is observed. It is possible that the slow-down 
during the kinematic dynamo regime is a consequence of a quasi-two-dimensionalization 
of the flow by the high rotation rate; this remains to be investigated fully. 

Figure [H] shows visualizations of energy densities and field lines at late times in 
Runs D5 and D7. The development of anisotropies in the velocity and magnetic fields 
can be observed in run D7, which has the highest rotation rate attained of f2 = 16. 
Indeed, as fl is increased the velocity field shows a tendency to develop columns, with 
mechanical energy concentrated in cylindrical structures aligned along O and with a 
larger component of than of f^. The velocity field in these column is helical, although 
in general the total mechanical helicity of the flow fluctuates in time around zero. These 
structures are observed before the magnetic field is introduced (although they persist 
as the magnetic energy grows) and seem to be the result of the Taylor-Proudman effect 
(see e.g. [3^ [35]). It is a trend observed through runs D5 to DIO (see Figure [15] for an 
example). 

Run D9 experiments with lowering the angle between the forcing's axis of symmetry 
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Figure 16. (a) Trace of the dipolc moment on the surface of the unit sphere (above) 
and ampHtude of the dipole moment as a function of time (below) for Run Dll. (b) 
Same quantities for run D12. 

and the axis of rotation to a = 20°, and the dipole becomes more difficult to excite (as 
evidenced by a smaller growth rate in the kinematic dynamo regime). Indeed, dynamo 
action could not be excited below 15° for the values of Reynolds numbers explored, as 
the resulting driven flow approaches axisymmetry. 

Run DIO is actually part of a set of experiments (Runs DIO to D12; again, 
see Table [1]) with forcing functions that are non-axisymmetric, and which may also 
inject mechanical helicity (Runs Dll and D12). Run DIO, having no net mechanical 
helicity, shows no big differences in the evolution of global quantities from the previously 
discussed runs. The dipole develops but its orientation wanders randomly, with some 
preferred orientation perpendicular to fl. In this case, axisymmetry is broken by the 
forcing directly instead of by a non-zero angle between an axisymmetry forcing and the 
axis of rotation. 

Runs Dll and D12 have non-axisymmetric forcing that injects mechanical helicity. 
The forcing for these runs is given by equation ( 1201) with coefficients 

■^3,3,0 = 5'C-3,3,0 = -^05 'C3,3,0<m<3 = 5,^-3,3,0<m<3 = -^o(l + i) , (21) 

with Fq = 1.7. In the presence of net helicity, dynamo excitation suddenly becomes 
much easier (as evidenced by a much larger growth rate of magnetic energy during the 
kinematic regime), and the ultimate saturation occurs at Em/Ev ~ 2: more magnetic 
than kinetic, with magnetic helicity, having sign opposite that of the mechanical helicity 
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Figure 17. Kinetic energy density and velocity field lines (left), and magnetic energy 
density and magnetic field lines (right) at the same time in the saturated state of run 
D12. The axes are aligned as in Figure [TBI 



inversely cascading to the large scales. As a result, the system is dominated by a helical 
magnetic field at the largest available scale. 

An interesting qualitative argument from mean field theory [36l [37] (which assumes 
large scale separations and often some form of periodic boundary conditions, as do all 
"a-effect" calculations) can be seen to anticipate this result as follows. From mean field 
theory , the induction equation for the mean magnetic field B (assuming there is no 
mean fiow U) is 



Here, a is proportional to minus the kinetic helicity of the fiow [361 EH ET], and /3 is a 
turbulent magnetic diffusivity. Dotting equation ([221) with the mean vector potential A 
(such as B = V X A) and integrating over volume, an equation for the evolution of the 
mean magnetic helicity Hm is obtained. 



where Em is the mean magnetic energy and Hj is the mean current helicity. As a 
result, if magnetic diffusion is neglected, the dynamo process injects into the mean 
(large) scales magnetic helicity of opposite sign than the kinetic helicity, and in the 
small scales magnetic helicity of the same sign. This effect has been observed before 
in numerical dynamo simulations with periodic boundary conditions [281 El]- The 
large scale magnetic helicity then inverse- cascades to the largest available scale in the 
system, while the small scale magnetic helicity is transferred to smaller scales where it 
is dissipated [39]. As a result, at late times the system is dominated by a large scale 
magnetic field with magnetic helicity of opposite sign to that of the kinetic helicity 
injected by the forcing. 

As seen in Figure [T6l the dipolar orientation in Dll seems to have a preference for 
being perpendicular to fl. For run D12, whose forcing function differs in its orientation 
to Q by 90°, the dipole orientation seems to remain in a single hemisphere (see Figure 



— = aV X B + /?V^B. 



(22) 



dHM 
dt 



(23) 
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[TBI) and its tip precesses about Q. Figure [T7| shows energy density and field lines in 
the saturated steady state of run D12. The axes are aligned as in Figure [161 As is 
often the case with helical flows, the geometry of the flow is more complex than in the 
non-helical runs. These runs in which a net sign of mechanical helicity is sustained 
by the mechanical forcing were continued for several thousands turnover times, and no 
reversals of the dipole moment were observed. The dipole moment seems to fluctuate 
around a preferred orientation with only short excursions of the dipole to the opposite 
hemisphere. 

5. Discussion and future directions 

By solving the mechanically-forced MHD equations inside a rotating conducting 
spherical boundary, we have found that a bewilderingly wide variety of dynamo behavior 
is possible for a magnetic Prandtl number of unity. The behavior is sensitive to 
mechanical and magnetic Reynolds numbers, to rotation rate, and more indirectly to 
the Ekman number. It is also sensitive to the geometry and strength of the forcing 
functions and their relation to the axis of rotation. We have only begun to explore this 
multidimensional parameter space. Note we have not made much effort to tailor the 
forcing functions we have chosen to model what the mechanical flows in planetary cores 
or stellar convective regions might be. Rather, we have been exploring dynamo behavior 
in the abstract, and feel somewhat overwhelmed by the variety of dynamo behavior that 
has been found. In this light, our attempt should be consider as an extension of dynamo 
simulations in periodic boundary conditions [271 EHl [311 SHI [29l [30] to consider the effect 
of boundaries and rotation, and not be compared with explorations of the space of 
parameters in realistic geodynamo simulations [lOl [HI [131 El [26] . 

What has become clear is that the wholly spectral methods we are using, while 
accurate, do not scale well into the parameter regimes of planetary and astrophysical 
dynamos, which involve many orders of magnitude between the largest length scales 
in the flows and the Ekman or dissipation scales that also play a role in the process. 
This limitation afflicts all numerical attempts to explore planetary and stellar dynamos, 
particularly in view of the low magnetic Prandtl numbers that are expected to prevail 
there, in simulations [291 [30], and in liquid-sodium experiments [HI [121 |l3l [HI US] . But 
the rapid multiplication of the number of terms to compute in the convolution sums 
in equations ([TOj) and (ITTl) provides a rather intractable limitation on efforts to use our 
code without modification at higher Reynolds numbers than the few hundred we have 
explored here. What seems to be called for is an exploration of the possibilities of 
using fast transforms (for spherical harmonics [46l [47] and possibly for spherical Bessel 
functions [^ ll9])to turn the code into a pseudospectral one in which the nonlinear 
terms are computed in configuration space rather than spectral space, as is commonly 
done for rectangular periodic boundary conditions [IT] which would increase available 
resolution by many orders of magnitude. Our future investigations will explore this 
possibility. 
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We also intend to replace the conducting boundary by an insulating but 
mechanically-impenetrable one that will permit protrusion of the magnetic field into 
the vacuum region surrounding the shell (see e.g. [2]). It may be that forcing the 
magnetic field lines to return to the interior each time they get near the shell that is 
now playing a dynamical role in what we are seeing would be different in the case of an 
insulating shell. This raises some conceptual difficulties, since the problem of matching 
MHD fields on to vacuum electromagnetic ones has not been thoroughly solved. For 
example, one approximation that has been used has been to match magnetic fields at 
a spherical surface to magnetostatic ones outside, which involve a magnetic field that 
is derivable from a scalar potential and for which there is no electric field. However, 
Maxwell's equations tell us that the tangential electric field must be continuous at any 
interface, and there is no reason why this tangential electric field should vanish or even 
be "small" immediately inside an insulating boundary of a conducting magnetofiuid. 
The magnetostatic approximation may be the best we can do, but it would be desirable 
to have more justification for it than we presently have. 

Finally, we need to devote attention to the geometry and strength of the forcing 
functions that are being employed and to study the effect of different forcing functions 
in dynamo action. Mechanical processes, convective and otherwise, are believed to 
power the dynamo in cases of geophysical or astrophysical interest. In the simulations 
discussed, we have no explanation, for example, why the "columns" or columnar vortices 
aligned along O form in the velocity field in runs toward the end of Tabled! It is clear the 
physical effect that is triggering their formation is the rotation alone, perhaps through 
a two-dimensionalization of the flow via the Taylor-Proudman effect [Ml [35| [8]. This 
is a different process than the conventional explanation for the formation of columns 
in planetary interiors, which involves both thermal convection and rotation [50], since 
thermal convection is completely absent in our incompressible MHD formulation. 
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